arXiv:1502.07944v2 [cond-mat.quant-gas] 10 Sep 2015 


Density-dependent synthetic magnetism for ultracold atoms in optical lattices 


Sebastian Greschner, 1 Daniel Huerga, 2,3 Gaoyong Sun, 1 Dario Poletti, 4,5 and Luis Santos 1 

1 Institut fur Theoretische Physik, Leibniz Universitat Hannover, Appelstr. 2, DE-30167 Hannover, Germany 
2 Instituto de Estructura de la Materia, CSIC, Serrano 123, E-28006 Madrid, Spain 
3 Institut fur Theoretische Physik III, Universitat Stuttgart, Pfaffenwaidring 57, 70550 Stuttgart, Germany 
*Singapore University of Technology and Design, 8 Somapah Road 487372, Singapore 
5 MajuLab, CNRS-UNS-NUS-NTU International Joint Research Unit, UMI 3654, Singapore 

Raman-assisted hopping can allow for the creation of density-dependent synthetic magnetism for 
cold neutral gases in optical lattices. We show that the density-dependent fields lead to a non-trivial 
interplay between density modulations and chirality. This interplay results in a rich physics for atoms 
in two-leg ladders, characterized by a density-driven Meissner- to vortex-superfluid transition, and 
a non-trivial dependence of the density imbalance between the legs. Density-dependent fields also 
lead to intriguing physics in square lattices. In particular, it leads to a density-driven transition 
between a non-chiral and a chiral superfluid, both characterized by non-trivial charge density-wave 
amplitude. We finally show how the physics due to the density-dependent fields may be easily 
probed in experiments by monitoring the expansion of doublons and holes in a Mott insulator, 
which presents a remarkable dependence on quantum fluctuations. 
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I. INTRODUCTION 

Orbital magnetism is crucial in condensed-matter 
physics. In particular, it plays a fundamental role in 
the integer and fractional quantum Hall physics m , as 
well as in related topics such as anyons [4] and topolog¬ 
ical insulators G3E]. Ultracold gases offer extraordinary 
possibilities for the controlled experimental simulation of 
quantum many-body systems 0. However, experiments 
are typically performed with neutral particles preclud¬ 
ing the direct quantum simulation of orbital magnetism. 
Interestingly, synthetic magnetism can be engineered in 
cold neutral gases, currently constituting a major topic 
in cold-atom research la Proper laser arrangements 
have been used to induce artificial magnetic fields and 
spin-orbit coupling mm In addition, recent experi¬ 
ments have demonstrated the creation of synthetic mag¬ 
netism in 2D optical lattices Ha na, and realized the 
Meissner-superfluid (MSF) to vortex-superfluid (VSF) 
transition |14j with bosons in square optical ladders [US¬ 
UI]. 

However, in all these experiments the fields created 
were static since there was no feedback of the atoms on 
the field. Nevertheless, such a dynamical feedback plays 
an important role in various areas of physics, ranging 
from condensed-matter m to quantum chromodynam¬ 
ics m ■ In order to experimentally simulate such dy¬ 
namical gauge fields in optical lattices various ideas have 
been recently suggested [20H28] . 

While in those kind of dynamical gauge fields it is cru¬ 
cial to discriminate particle and gauge fields degrees of 
freedom, gauge fields determined by the particle den¬ 
sity appear in a variety of problems in condensed-matter 
physics, including effective field theories for the fractional 
quantum Hall effect [29] and their extension to general¬ 
ized statistics in one dimension m, spin models of quan¬ 
tum magnetism [31 , and chiral solitons [32] . Indeed, the 


atomic back-action on the synthetic gauge field in cold 
gases experiments is expected to lead to these last type 
of excitations in Bose-Einstein condensates [33] . 

From the experimental standpoint, Raman-assisted 
hopping of cold atoms [12j Ena n3nzi can be used 
in ID optical lattices to realize occupation-dependent 
gauge fields that result in effective anyon statistics and 
thus a clean experimental realization of a ID anyon- 
Hubbard model [sum. This model presents a wealth of 
new physics, including statistically-induced phase transi¬ 
tions [34], novel superfluid phases [35] . asymmetric mo¬ 
mentum distributions [36], and intriguing dynamics [33- 
[39]. 

In this paper we show how a Raman-laser scheme can 
be employed for the realization of density-dependent syn¬ 
thetic magnetism (DDSM) in one and two dimensions 
and study its effects in ladder and square lattice geome¬ 
tries. In particular, we are interested in a system de¬ 
scribed by the following Hamiltonian 

n = -Z [tA +ex e^ n 'b r + t y bl +e b r + H.c. 

r 

+ — y ] n r (n r — 1) — (i y ^ n r , (1) 

r r 

where b r (bj) is a bosonic annihilation (creation) opera¬ 
tor acting on site r = (i, j) of the lattice, and n r = b\b r is 
the number operator. As it will be discussed below, the 
experimental implementation leads to a three-body hard¬ 
core constraint on the onsite occupation, i.e. n r = 0,1, 2. 
The first term in 0 accounts for the hopping of bosons 
along the two directions of the lattice, defined by lattice 
vectors e x = (1,0) and e y = (0,1), while the second and 
third terms account for the usual onsite Hubbard inter¬ 
action and the chemical potential which fixes the total 
density of the system, respectively (see Fig. [l] (a) and 
(b)). As it will be shown below, the density dependent 
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FIG. 1: (a) Sketch of the density-dependent Peierls phases of the model 0 on a ladder (see also Eq. ( p~2| )) and (b) of the 
2D-square lattice, (c) Scheme of the creation of a density-dependent Peierls phase using Raman assisted hoppings, (d) Raman 
assisted hops (i)-(iv) as discussed in the main text. 


Peierls phase of the hopping amplitude (e 1(?!)rnr ) can be 
chosen in such a way that an effective net-magnetic flux 
per unit-cell is created. In particular, we will concentrate 
in the case were the phases depend only on the position 
in the ^/-direction, i.e. p r = pj. Due to the operator 
nature of this phase, quantum fluctuations of the density 
will crucially affect the effective magnetic flux. In this 
work, we demonstrate that DDSM has important con¬ 
sequences for bosons in two-leg ladders and 2D square 
lattices, leading to a non-trivial interplay between chiral¬ 
ity and density modulations. 

The structure of the paper is as follows. In Sec. [Tl 
we comment on the realization of DDSM using Raman- 
assisted hopping. In Sec. m we analyze the conse¬ 


quences of DDSM in optical ladders, whereas in Sec. IV 
we focus on the case of 2D square lattices. In lad¬ 
ders, this interplay results in a density-driven Meissner- 
superfluid (MSF) - Vortex-superfluid (VSF) transition 
with a non-trivial density imbalance between the legs. 
In square lattices DDSM induces a similar transition be¬ 
tween a non-chiral superfluid (SF) and a chiral super¬ 
fluid (CSF), both presenting a non-trivial density-wave 
amplitude. Section[V]is devoted to the dynamics of parti¬ 
cles and holes, which are crucially affected by the DDSM, 
as illustrated by the expansion of doublons and holes in 
a Mott insulator (MI), which presents an intriguing de¬ 


pendence on quantum fluctuations. Finally in Sec. VI 
summarize our results. 


II. DENSITY-DEPENDENT FIELDS 

In this section we propose a possible experimental 
scheme for the realization of DDSM. First, we briefly re¬ 
view the proposal for the creation of a density-dependent 
Peierls phase in one dimensional lattices, as described in 
Ref. [35] . which is the key ingredient for the realization 
of DDSM. In following subsections, we discuss how this 
scheme naturally extends to higher dimensional lattices 


and how it may be adjusted to effectively reproduce the 
density-dependent Peierls phases of Model 0. 


A. Two-component system 

We consider a bosonic species with two internal states, 

| A) and | B). As shown below for the specific case of 
87 Rb, we may choose | A) = \F — l,mp = —1) and 
| B) m |_F = 2, m y = —2). A detailed discussion of other 
species can be found in the supplemental information of 
Ref. [35] . Both components are confined to the lowest 
band of a tilted ID optical lattice along the x axis, of 
spacing D, and depth Vo = sEr, with Er = h 2 ir 2 /2m D 2 
the recoil energy. The Hilbert-space of a single lattice 
site thus constitutes of empty sites (0), single occupied 
sites (A) or (B), doubly occupied sites (AA), (BB) or 
(AB), etc. Without tilting there is a hopping rate J to 
nearest neighbors. The lattice tilting induces an energy 
shift A from site to site as shown in Fig. 0(c). 

We denote as w(x—jD ) the Wannier function at site j. 
Due to the tilting, it is convenient to use Wannier-Stark 
states. For J « A, the Wannier-Stark state centered 
at site j may be approximated as ^ j{x ) — w(x — jD ) + 
^[w(x — (j + 1 )D) — w(x — (j — 1 )D)] [40]. The 3D on¬ 
site wavefunction at site j is 4>j(r) = !pj(x)(p(y, z\ where 
c p(y,z ) is given by the strong transversal confinement. 
For simplicity we assume below tp(y,z) ~ w(y)w(z). 

On-site interactions between atoms in states a and 
/3 (for a, /3 = A, B) are characterized by the cou¬ 
pling constant U a $ = 47r ^°^ f d 3 r\w(r)\ 4 , with a a ^ 
the corresponding scattering length. For a sufficiently 
deep lattice, the evaluation of the on-site interactions is 
simplified by means of the harmonic approximation [7k 
<F(r) (y / 7r/) _3 / 2 e _r2 / /2 , where l = Ds _1 / 4 /7r. Using 

this approximation we obtain: U a/27t 5 / 2 ,s 3 / 4 . 

As shown in Ref. [35] the scheme may as well be realized 
with fermionic species; then inter-species on-site interac- 
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tion Uab is possible. 


B. Raman-assisted hopping 


No direct hopping occurs since J< A, U a ^. Raman- 
assisted hopping is realized with the set-up of Fig. [l] (c) 
formed by four lasers Lj— with Rabi frequencies 
Qj = ^ wave vectors kj, and frequencies ujj. L\ ,4 

have linear polarization and 1/2,3 circular cr_ polarization 
and couple states | A) and | B) far from resonance. \B) 
is just affected by lasers L\ 4 due to selection rules. Al¬ 
though both Z/2,3 and L\ ,4 couple to | A), the coupling 
with Li ? 4 can be made much smaller than that of 1 / 2,3 
(for a detailed discussion see the supplemental informa¬ 
tion of Ref. [35]). Hence we may assume below that | A) 
is just affected by L 2 , 3 . 

Following Ref. m , we evaluate the Raman-assisted 
hopping, J nm , given by lasers L n=lj2 and I/ m=3?4 , from 
site j to site j - hi: 

Jnm = J d 3 r$ j+1 (r)*e MkBm - r $ j (r), (2) 

where = <p n - (f> m , Sk nm = k„ - k m , and = 
| ^ (i ar g e ) detuning to the one-photon 

transitions. For J<A and s 1, we may approximate: 

J nm c i J sin .. 

(3) 

Note that 7 ^ 0 is necessary to establish a significant 
assisted hopping mmm- Each laser pair couples a 
different Raman transition (see Fig. 0 (d)): 

(i) J 23 characterizes the hopping (A, 0 )—)>( 0 ,A), which 
is accompanied by an energy shift A E = — A. We 
hence demand u) 2—^3 = —A and the transition 
amplitude is given by V23 — \ including the 
appropriate Clebsch-Gordan coefficients for the spe¬ 
cific case of 87 Rb. 


(ii) (A,A)—>>((), AB) is given by J 24 , being characterized 
by A E = —A + Uab ; we impose co >2 — ^4 = —A + 
Bab + B , with U <C Bab, A and the amplitude 


U 2 4 = 


1 O2O4 

76 ^“ 


(iii) J 13 is linked to the hop (AB, 0 )—»(A,A); the energy 
shift is A E = —A — Bab ; we demand uji — ujs — 
-S-Uab-U. Vi3 = 


(iv) (AB,A)—»(A,AB) is given by J 14 ; the energy shift 
is A E = —A; we impose uq — = —A. The 

transition amplitude is given by Vi 4 = | _ 


The frequencies ujj are chosen such that they compen¬ 
sate the lattice tilting, and hence no Bloch oscillation 
is induced in the rotating frame. In this frame process 
(ii) is accompanied by an energy shift [/, (iii) by a shift 


—17, and (i) and (iv) have no associated shift. U may 
be hence understood as an effective on-site interaction 
energy. Alternatively, these energy shifts are compatible 
with an on-site interaction Bab, and an effective nearest- 
neighbor interaction V = Bab — U. We will return to 
this point when discussing the extension to 2D lattices. 

Note that processes (i) and (iv) are energetically degen¬ 
erated, but they may be addressed with different lasers 
due to selection rules. This point constitutes the major 
drawback of the proposal of Ref. [34]. In that proposal, 
a single component, A, was considered, and process (iv) 
was of the form (AA,A)—^(A,AA), which cannot be re¬ 
solved from process (i). As a result, in the scheme of 
Ref. [34], both the combination of L 2 and I/ 3 , and of L\ 
and I /4 address both (i) and (iv), preventing the realiza¬ 
tion of the desired density-dependent Peierls phase. The 
two processes may be just discerned by considering a very 
small detuning S < Baa, A, which would be accompanied 
by very large spontaneous emission losses. 


C. Spurious processes 

Undesired spurious processes are in principle possible: 

(v) (A,0) -> (0,B); AE = -A 

(vi) (A,A) -► (0,AA): AE = -A + Uaa 

(vh) (AA,0) -► (A,A): AE = -A - U A a 

(vih) (AB,A) (B,AA): AE = -A + SU, with SU = 
(Baa ~ Bab ) 

(ix) (AA,B) -► (A,AB); AE = -A - SU 

Process (v) is just possible with J 24 or J 13 . But these 
laser combinations are (quasi-)resonant with — A±Uab- 
For Bab W, with W the width of the Raman reso¬ 
nance (typically of the order of 50 Hz [40]), process (v) 
is far from resonance with either J 24 or J 13 . To neglect 
the (vi) and (vii) processes one needs Baa W. In con¬ 
trast, to avoid (viii) and (ix) one must demand SU W . 
The latter condition is certainly more strict, but may be 
attained in experiments as shown in the supplemental 
information of Ref. [35j. 


D. Effective ID Hamiltonian 


W r n -r, lirnr |Oi||0 4 | _ |0 2 ||0 3 | _ \Qi\\Q 3 \ _ I ^ 2 11^41 _ 

We assume — - 3 - ^ - 2v ^ - 

O 2 , Oi = |f 2 i|e _ q and 2,3,4 = O j and obtain the 

/~\ 2 /— 2 

transition amplitudes U 23 — ^~, V 24 = v2y, V 13 = 
and Vi 4 = 2^-e -1 ^. Note that an additional 
factor y/2 is used to mimic bosonic enhancement. We 
denote as Cj the bosonic operator corresponding to the 
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Fock-state manifold {|0), |1) = \A), |2) = \AB)}. As¬ 
suming ki ? 2 = ke y , and = ke x , and kD = 7r, then 


'H 




cje'^cj+i + H.c. 


U 


7E n i( n i _1 ) 


(4) 

with rij = cjcj, and t = ( ) T Typical values of the 

Raman-assisted hopping rate, £, are of the order of few 
tens of Hz 00]. Note that the factor (—l)- 7 , which results 
from the x projection of $k, may be easily eliminated by 
redefining the bosonic operators in the form: b^i = 
^ 4 Z+i = Qz+i, ^ 4 Z +2 = —C 4 /+ 2 , ^ 4 Z +3 = ~c 4 z+ 3 , with l an 
integer. In this way we obtain the ID model: 


n = -tJ2 [b]e^b j+ 1 +H.c. 

3 


^XA( n J -1 ) ( 5 ) 

3 


E. Density-dependent gauge fields in 2D lattices 

For a 2D square lattice or ladder, one may proceed as 
in Refs. 02 nu, assuming assisted hopping along x, and 
natural hopping along y. This is however problematic, as 
one can clearly understand from the alternative picture 
mentioned above (section II B), in which the on-site in¬ 
teractions remain characterized by Uab , but an effective 
nearest-neighbor interaction, V, is induced along x. In 
contrast, along y there is no nearest-neighbor interaction. 
Although this asymmetric extended-Hubbard model may 
have interest in itself, it is not the model to be explored 
in this work. 

An effective model with only on-site interactions and a 
density-dependent gauge demands both directions to be 
Raman-assisted. Following the same arguments as above, 
we evaluate the assisted hopping given by lasers n and m 
from a site r = (D x r x , D y r y ) to the site r + with 

ej the unit vector along the j = x,y direction, and Dj 
the lattice spacing along that direction: 

4ft = (^) e 1 *™ sin 

( 6 ) 

where A j and Jj are, respectively, the tilting and the 
hopping without tilting along the j direction. 


F. Four-laser arrangement 

We first consider the same arrangement of four Raman- 
lasers as discussed above. We assume A x = A y = A, 
and J x = J y = J. For ks = k 4 = jy-e x , k 2 = and 

k i = 

J^ x ’ v) (r) = (y^K~) e^^ + r v+ i,2) x uy (7) 


with = 1, and A^ 14 = A^ 14 = 

g^Gb+V 2 ) cos(0/2). Assuming Q 2 ,3 = £2, ^ 1,4 = V2Q, 
4 = 0? we obtain the effective Hamiltonian: 


n = 


+ 


-fyy- 1 ) r * +r » ^2 [ct +e Jj(r y ,n r )c T + H.c. 

r j=x,y 

^y]n r (n r -1). (8) 

r 


with f x (r y ,n r ) = e ^ nr , and f y {r y ,n r ) = 

cos(0n r )e 1 ^ r ^ +1 / 2 ^ )nr . As for the ID problem, the 
factors (—1 y*+iy originate from the projection of 
all Sk nm along x and y in order to achieve as¬ 
sisted hopping along both directions. They may be 
readily eliminated by introducing the transforma¬ 
tion b 4 Ux y Sx An y -\-Sy = Sy)c4 nx + Sx An y +Sy 5 where 

c(s x , s y ) = —1 if ( 83 , + s y ) mod 4 > 1 and e(s x , s y ) = 1 
otherwise. With this transformation we obtain: 

K = ~tJ2 J2 [t>t+ ej fj{r y ,n r )b T + K.c. 

r j=x,y 

+ ^yy r (n r -i). (9) 


With a four laser arrangement it is hence possible to 
create a density-dependent Peierls phase in the 2D lat¬ 
tice. However, since S k nm must project on both x and y 
directions, there is an additional dependence of the hop¬ 
ping modulus along y on cos(0n r ). Similar dependences 
of the hopping modulus appear for other choices of ki... 4 . 


G. Six-laser arrangement 


A model in which the tunneling modulus does not 
depend on the occupation may be attained by adding 
two additional lasers. This six-laser arrangement al¬ 
lows as well for a more flexible realization of density- 
dependent Peierls phases. We consider A x ^ A y , and 
two additional lasers I/ 5 ?6 , with cj 5 = uj\ T A y — A^, and 
ujq = co >2 + A y — Aa,. In this way, the hopping processes 
(i) to (iv) along y are produced, respectively, by Jq%, 

J53, and J54. If \A X — A y \ IF, the y hops may be hence 
addressed independently from those along x. We impose 
Sk^Dx = 7 r for n = 1,2 and m = 3,4, and Sky ,rn D y = n 
for n' = 5,6 and m = 3,4. A possible example is given 
by k 3)4 = -J^e x , k 2 ,5,6 = and ki = ^ e y . We 

choose £2 2 /\/3 = r^i/2 = I)', 12 3 /\/3 = I2 4 /2 = 12 and 
n 6 /V3 = ^ 5/2 = £ 2 " and after eliminating the factors 
(_1 y x +j y as above, we obtain the Hamiltonian 0 with 
tx = ai ‘b ty = ail( -l 4> y = 4>r y - Although 

also the system 0 exhibits DDSM, in the following we 
will for simplicity focus our analysis on Hamiltonian 0 
with pure density-dependent Peierls phases. 
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FIG. 2: (a) and (b) Phase diagram for a ladder with t y = t x , 0i = 0 .87T, and 02 = 0 as function of t x /U and chemical potential p 
computed with the density matrix renormalization group (DMRG). The color code indicates (a) the particle density imbalance 
between the legs An and (b) the chiral current J c (obtained from simulations with L — 24 rungs). Solid lines mark the MI 
with p — 1, and (very narrow) with p — 1/2 and 3/2 (extrapolated to the thermodynamic limit from systems with up to L = 96 
rungs). The dashed line denotes the MSF-VSF transition, (c) An (dashed) and J c (solid) for the same parameters as in (a) 
and (b) but U — 0 and L = 48. (dl-5) Typical particle density and current configurations for U — 0 and (dl) p — 0.1, (d2) 
p — 0.63, (d3) p — 1.25, (d4) p — 1.46, (d5) p — 1.77. The size of the circles is proportional to the onsite-density, he arrows 
encode the strength of the local currents. 


III. DDSM IN LADDERS 


Ladders with static fields (i.e. Hamiltonian 0 with 
density-independent Peierls phases) have been recently 
realized in several experimental groups mm and stud¬ 
ied theoretically as well mEMi]. 

An important observable in this context of density- 
independent fields, but also for DDSM, are currents m • 
From the continuity equation 



([H,n r ]) 


( 10 ) 

<S> 


we can define the current J(r —> s) from a site r to a 
neighboring site s. 

At low fluxes the system is a MSF [14], characterized 
by the absence of rung currents. At a critical flux, which 
depends on t x /t y and interactions, currents penetrate the 
rungs, form vortices, and the system becomes a VSF. 
Figures of typical current configurations of the VSF and 
MSF phases may be found in Fig. [2] (d) anticipating the 
discussion below. The MSF-VSF transition is signaled 
by a cusp in the chiral current J7/ Static = J _ jF 2 statlc 5 
with 

^static = }_ ^ (6 t , e -i* 6j . +li . _ H.C.) (11) 

3 


(in units of the hopping t x ) the leg currents of the ladder 
with density independent static magnetic Peierls phases. 
For bosonic systems with a finite interaction V, vortices 
may form crystals of a commensurate vortex-density py, 
which is a not-conserved quantity measuring the number 
of vortices per system size. Such so called vortex-lattice 
(VL pv ) phases have been studied in weak coupling 
regime m as well as in numerical calculations [46]. 


In the following we study the DDSM in a ladder geom¬ 
etry shown in Fig. [l] (a). In this situation given by the 
following Hamiltonian 


Uc = — tx V [ b m,i ei4,inj ' lb i,i +H.C.' 

3 

-^E[Vl,2^ 2 " 3 ’ 2 ^2+H.C.' 

3 


-ty E [kj,2^7,1 + Hc - 
i 



( 12 ) 


A. MSF and VSF phases with imbalanced density 


In the limit of strong on site repulsion U J 
model (Jl2| with density dependent phases may easily be 
mapped onto a system of hardcore bosons without a flux 
for fillings 0 < p < 1. For 1 < p < 2 we may con¬ 
sider doublons \2)i j on top of a uniform Mi-background 
n,j uk j as hardcore particles, which, however, now ex¬ 
perience a finite flux 0 = 0i — 02 , such that the effective 
Hamiltonian in this limit may be written as 


hardcore = E [ C U/^ + H ' C 

hJ 

[ c ],2 c id + H* c * , 

3 


(13) 


with Cij ( 4 .) being the creation (annihilation) operator 
of a doublon on site (z, j). Thus in the strongly interact¬ 
ing regime U J, model © is expected to reproduce 
the physics of hardcore bosons in a magnetic static field 
exhibiting MSF and VSF phases as discussed in detail 
in [45] . 








































6 


(a) (b) 




FIG. 3: (a) Phase diagram for a ladder with U = 0, 0i = 7r, and 02 = 0 as function of t y and p (in units of t x = 1). As 

discussed in the main text for large densities a VSF is realized, while for low densities the system is in a MSF phase, separated 
by a commensurate-to-incommensurate (with respect to the vortex density pv) phase transition (dashed line). Dotted lines 
denote lines of constant particle density p — 1/2,1,3/2 and 7/4, while solid lines indicate the gapped phases, (b) Phase diagram 
as function of p/U and 0i = — 02 = 0 for U — t x and t y = 10 t x computed with DMRG. The shaded areas denote gapped phases 
of p — 1/2,1, 3/2 and 2 filling - the white area correspond to MSF phases. As a characteristic feature of the density-dependent 
fields one observes a sequence of direct transitions between the gapped phases resulting in a macroscopic jump of density 
around 0 = 7r, 2tt/3 and 7t/2 (see text). All gapped phases exhibit finite rung-string order Orso and rung-parity order Orpo 
as defined in the main text. The inset shows Orpo (solid line) and Orso (dashed line) for p — 1/2. 


A qualitative insight on the physics induced by the 
occupation-dependent Peierls phases away from that 
limit is provided by a simple mean-field decoupling of 
the tunneling terms (between neighboring sites r and 
r') in Eq. (Jl2|: Since (6j) 3 = 0, then b\.,e l ^ nr b Y = 
&J,( 1 + (e^ — 1 )n r )b r . Using the decoupling b\,n Y b r ~ 
2ft (r', r)(n r — n Y ) + 2 n Y b\,,b Y we obtain 

b\,e 1(t>n *b r + H.c. ~ [(1 + 2n r (e i0 - 1 )6+ / 6 r + H.c.] (14) 

+ [2ft(r / ,r)(e 1( ^ > - 1) + c.c.](n r - n r ), 


with n Y = (n r ), and ft(r',r) = (b l,b r ). The first term at 
the right hand side (rhs) of Eq. ( |l4| ) results in an effec¬ 
tive Peierls phase. A density-dependent effective flux is 
hence given by the phase accumulated when encircling a 
plaquette. For model (12) the effective flux is uniform 
for homogeneous n^j = f0. The second term at the rhs 
introduces a shift of the local chemical potential, which 
is leg-dependent for 02 ^ — 0i in Eq. (12). We hence 
expect from this simple argument that the occupation- 
dependent Peierls phase introduces an interplay between 
density-dependent fields and density imbalance between 
the legs. We employ below density matrix renormaliza¬ 
tion group (DMRG) |47 calculations to confirm this in¬ 
sight. In these calculations we use system sizes up to 100 
rungs with open boundary conditions keeping up to 1000 
matrix states. 

As in Ref. m we monitor the chiral current, J c = 
Ji — J 2 , where the leg currents are now defined as 


Ji = { , • u - H.c.) (15) 

3 


in units of t x /h, with TV/ the number of bonds along 
the leg. The density-dependent effective flux induces a 
characteristic dependence of J c on chemical potential, 
/i, as shown in Figs. 2] (b) and (c) for 01 = 0.87T, and 
02 = 0 (similar results occur for other parameter val¬ 
ues). In addition to MI phases, we observe two different 
SF regimes. For a given t x /U there is a critical /i, at 
which J c presents a cusp (see Fig. [2jc)), indicating a 
MSF-VSF transition induced by the increasing effective 
flux for growing lattice filling p. 

This transition is as well characterized by a kink in the 
equation of state p(p) which signals a change in the num¬ 
ber of gapless modes of the system: The MSF phase has 
a gap in the antisymmetric (or naively “Vortex”-) sec¬ 
tor m, while the symmetric (“charge”) sector remains 
gapless. The number of gapless modes is also reflected 
by the central charge c, which is c = 1 in the MSF phase 
and c = 2 in the VSF phase. We verify this by numer¬ 
ically extracting c from the scaling of the entanglement 
entropy [45] . 

Figures |2ja) and (c) also show that the occupation- 
dependent Peierls phase leads to a marked density im¬ 
balance, An = 2(n 2 — ^i)/(^2 + fii). Three important 
points should be noted. First, although \An\ is particu¬ 
larly large in the MSF, it is non-vanishing as well within 
the VSF. Second, An presents a kink at the MSF-VSF 
transition. Third, although An results from the explic¬ 
itly broken symmetry between the legs in Eq. (fl2|, its 
sign depends non-trivially on p or p. Figure [2|cjshows 
that An may change its sign going through a balanced 
point, An = 0. 
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B. Strong rung-coupling limit 


As for the case of static magnetic fields where the 
MSF-VSF transition has been explored for a fixed flux 
in Ref. [15j as function of t y /t x , also for the DDSM the 
rung hopping strength constitutes an important degree 
of freedom. In Fig. [3] (a) we study the phase diagram, in 
particular the commensurate-to-incommensurate MSF- 
VSF transition, on t y for U = 0, pi = 7r, and p 2 = 0. 
Interestingly we basically observe two different regimes: 
For small interchain couplings t y /t x < 1 the MSF-VSF 
boundary (dashed line) is located close to unit filling, 
however, as t y /t x > 1 it shifts quickly to larger densities 

p rsj 3 / 2 . 

For p i—02 7 ^ tt above a critical value of t y /t x the 
VSF phase may vanish. In the strong-rung coupling 
limit t y /t x 1 several gapped band insulating phases at 
commensurate fillings p = 1 and 3/2 are stabilized. The 
extent of the Mi-phases strongly increases with t y /t x 
as discussed below. For the parameters of Fig. [3] (a) a 
Mi-phase at p = 1/2 is suppressed. Additionally one 
may observe a gapped charge density wave phase at 
filling p = 7/4 (see [45] for a detailed discussion of similar 
phases at 1/4 filling for density-independent synthetic 
magnetism). Apart from the SF-phases also the MI- 
phases may be of Meissner-MI (in Fig. [2|a) for p = 1/2 
and p = 1) and of vortex-MI (for p = 3/2) types as 
discussed in Ref. m- Both exhibit a mass gap, however, 
the vortex-MI still has gapless mode (i.e. the neutral 
gap in the manifold of constant particle number van¬ 
ishes) while the Meissner-MI phase is completely gapped. 


For t y U , t x the ladder reduces to an effective rung- 
chain model with intriguing novel features due to the 
density-dependent Peierls phases. We may then map 
to rung-states | N) with a fixed occupation, of N = 
0,1, 2, 3,4 particles, on each rung. For the particular 
case of pi = — p 2 m <P at t x = 0, the ground-states of the 
decoupled rungs are the rung states 


| 6 > 

|i> 

| 2 ) 



U 

5 ty p ^ t y H - vr * 


U 


U 


5 ty 4 * < P < , 


where the notation denotes the rung state with 

n 1 (77,2) particles in the upper (lower) leg. At low p (in 
the vicinity of p ^ —t y ), for which N > 2 are irrelevant, 
the effective rung-chain model becomes of the form: 

n = -t x J2 [b){ 1 - sin 2 (<V2 )Nj)B j+1 + H.c. 

3 

+ ( 16 ) 

3 3 



FIG. 4: (a) Density imbalance An and (b) chiral current 

J c function of pi and p 2 for t x — t y , U — t x and p = 1.25 
as obtained by DMRG calculations. In addition to the MSF 
and VSF phases a small vortex lattice phase at vortex-density 
pv = 1/2, VL 1 / 2 -SF may be observed. Dashed lines indicate 
the phase boundaries from the VSF to the MSF and VL 1 / 2 -SF 
phases. 


where Bj are bosonic operators in the space {|0), |1), 12)}, 
and Nj = B^Bj. Note that inter-rung hops |0) 11) 

11)10) and 11)11) o |0)12) have an amplitude t Xj whereas 
11)11) |2)|0) and |2)|1) o |1)|2) have an amplitude 

t x cos 2 ((p/2). The latter rate vanishes for p = n. As a 
result a direct transition occurs for finite t x between the 
gapped phases at fillings p = 1/2 and 1, i.e. N = 1 
and 2, with an infinite compressibility and a macroscopic 
density jump (Fig. [3^b)). Similarly, direct transitions 
occur between gapped phases with p = 1 and 3/2 (at 
p = 27t/ 3), and 3/2 and 2 (at p = tt/2). 


The presence of the density-dependent phases results 
in a broken space-inversion symmetry, since the ampli¬ 
tudes of 11)| 1) 12)10) and 11)11) |0)|2) are not equal. 
As discussed in Ref. m the broken space inversion sym¬ 
metry may result in the exotic situation of simultaneous 
presence of both nonlocal parity- and string-order in the 
insulating Mi-phases. The MI phase of a usual ID Bose 
Hubbard model is characterized by a finite hidden par¬ 
ity order due to bound particle hole pairs that has been 
observed in experiments with single site resolution [45] , 
A non-vanishing string order, but vanishing parity or¬ 
der, characterizes the Haldane insulator, predicted in 
polar lattice gases mm and bosons in frustrated lat¬ 
tices [52] . The explicit expressions in the effective rung- 
state model may be borrowed from the corresponding 
orders of a spin 5=1 chain m- We define the rung- 
parity-order 0 RPO = lim| i _ i |_ >00 ((-l) S *<'=<i SNk ) (with 
SNk = N — Nk) and the rung-string-order Orso = 
lini|j_j|_ ) . 00 ((5A^(— SNk SNj). For 0 = 0 the 
Mott phases of the rung-chain model (16) present fi¬ 
nite rung-parity-order Orpo but vanishing rung-string- 
order Orso • Dae to the density depended phases the 
Mott rung phases acquire a simultaneous finite Orpo 
and Orso , as may be seen in the inset of Fig. 1(b)- 



















C. Symmetries and vortex-lattice phases 


Density independent magnetic fields are up to a gauge 
transformation completely defined by the net flux per 
unit-cell of the lattice. Due to its operator-nature this 
is not true for the case the density dependent Peierls 
phases. Indeed as may be seen in Fig. [4] the phase di¬ 
agram may significantly depend on the values of both 
phases fa and fa of model (fT2l) . While the MSF-VSF 


phase boundary mainly just depends on the total effec¬ 
tive flux fa— fa, only in the vicinity of (j)\ ~ —fa ~ 7t/ 2, 
where also density imbalance An vanishes, we observe a 
vortex-lattice phase at vortex-density py = 1/2 (VL^ 
phase). Apart from the characteristic staggered pattern 
of the currents as shown for the case of static magnetic 
fields e.g. in [46], the VL]y 2 phase may be discriminated 
from the VSF phases by the calculation of the central- 
charge, which is c = 1 in vortex-lattice phases. For strong 
phase-imbalances fa = 0, fa = (j) as in Fig. [2] no vortex- 
lattice phases are observed. 

Note that the phase diagram is symmetric with respect 
to inversion of the phases I\ : ( fa, fa) — (—fa, — fa) and 
exchange of the two legs of the ladder / 2 : (fa, fa) 
(fafai) for fa, fa G [0,27r). The density imbalance An 
(Fig. 0] (a)) is (anti)symmetric with respect to ii(/ 2 ). 
The chiral current (Fig. [4] (b)) is an antisymmetric quan¬ 
tity under both A and J 2 . 


IV. TWO-DIMENSIONAL SQUARE LATTICES 

We now extend our study to the effect of DDSM to 
two dimensional square lattices. Equivalently to the lad¬ 
der case, we show that the occupation-dependent Peierls 
phase induces a non trivial interplay between the density- 
dependent phases and density modulations. As a first 
approach, we focus on the limiting case of 7r-phases, i.e. 
fa = jbr, for which the Peierls phase in Hamiltonian 0 
takes the simpler form, 

e i0 r n r = ( 17 ) 

Despite this simplification, the Peierls phase of the hop¬ 
ping along the x-direction still rends the 2D model 0 to 
be intrinsically frustrated and thus highly non-trivial to 
approach from a computational perspective. In the fol¬ 
lowing, we will use the composite boson mean-field the¬ 
ory (CBMFT) [53j [54], which is a useful tool to unveil 
strongly correlated phases of spin and boson lattice mod¬ 
els where other methods face significant problems. 

CBMFT is based on the use of clusters of the original 
degrees of freedom as the basic degrees of freedom that 
contain the necessary quantum correlations to describe 
the phases emerging in the system under study. In prac¬ 
tice, we tile the lattice into clusters of equal size, in such a 
way that each site r of the original 2D lattice belongs to a 
unique cluster. The tiling is performed preserving most 
of the symmetries of the model. Each quantum state 


of each cluster can be represented by the action of a cre¬ 
ation composite boson (CB) over a CB vacuum. Being the 
mapping relating the original bosons {frj, b r } to the new 
CBs canonical [53], one can rewrite (fTl) in terms of CBs 
and approach it by standard many-body techniques, with 
the advantage that short-range quantum correlations are 
exactly computed by construction. 

Here we will use the CB Gutzwiller ansatz, a simplest 
product of uncorrelated cluster wave functions, 


i$)=n Agi 0 * 75 ) = n i g ) R 


(18) 


R 


R 


where aj^ g is the creation CB associated to the clus- 

ter R in the state |g)R = Un ; |ii)r, where n refers 
to a cluster configuration in the occupation basis. The 
amplitudes Un ; are then determined upon variational 
minimization of the energy. In the homogeneus case, i.e. 

fR\ 

Un — Un, this variational determination is equivalent 
to exactly diagonalize a unique cluster with open bound¬ 
ary conditions and a set of self-consistently defined mean- 
fields acting on its borders [ 


The CB Gutzwiller ansatz (18) allows to compute ob¬ 


servables and order parameters in a systematic way. In 
particular, the energy obtained is variational, and the 
ground state phase diagram can be obtained by monitor¬ 
ing the ground state energy and its derivatives. In addi¬ 
tion, low lying excitations over the ground state can be 
analyzed within the CBMFT framework self-consistently 
[53] . Nevertheless, this analysis is out of the scope of the 
present work. 

We define a (0, tt) charge density wave (CDW) order 
parameter, pcdw = ^2 r (4>|n r \fa/N, which com¬ 

puted with an homogeneus CB Gutzwiller ansatz |g) 
takes the form, 


1 


Pcdw = 


E 


LxL y ren 


e (g|n r |g)- 


(19) 


and the bond-chiral order parameter, 


( 20 ) 


<r,r'> 


where is the number of bonds, and the currents J r y 
are defined through the continuity equation (10). 

Figures [5] (a) and (b) show CBMFT results of the 
ground-state phase diagram of Eq. 0 with 7r-phases and 
clusters of size L x xL y = 2x2, 4x2. These sizes preserve 
the periodicity imposed by the Peierls phase with effec¬ 
tive 7T-flux. In order to enhance the non-trivial hopping 
of bosons along the ^-direction we have set t x = 2 t y . 

As we can see in Figs. [5] (a) and (b), the system 
presents the usual MI lobes of integer density for small 
values of the hopping t x . For bigger values of the hop¬ 
ping, the ground state presents superfluid order, char¬ 
acterized by a nonvanishing condensate density, po = 
($|6j c _o^k=o|^)/A" (not shown). In particular, a SF 
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FIG. 5: (a) and (b): Ground state phase diagram of model 0 in 2D with 7 r-phases 0 r — jit and t x — 2 t y computed with 

CBMFT with clusters of size L x x L y — 2 x 2 (black lines) and 4x2 (dotted grey lines). Solid lines mark the boundaries of 
the MI, while the the onset of a finite chiral bond order parameter, 77 , signaling the CSF-SF transition, is marked with dashed 
lines and the dashed curve marks. The color code indicates (a) the charge density wave order parameter pcdw and (b) the 
chiral bond order parameter 77 , as defined in the main text. Note that the CBMFT results remain stable under increasing of 
the cluster size, (c) Chiral bond order parameter 77 (solid line) and charge density wave order parameter pcdw (dashed line) 
for a cut in the phase diagram along t x /U = 0.1. (d) Typical current and density configuration of the CSF phase. The size of 
the circles is proportional to the onsite-density, the lengths and widths of the arrows encode the strength of the local currents. 


phase with modulated density and vanishing bond-chiral 
order emerges for p < 1, while for p > 1 the ground state 
is a CSF characterized by having nonvanishing bond cur¬ 
rents in a pattern of fully stacked checkerboard pattern 
of vortices and anti-vortices (Fig. [5] (d)). Notice that the 
CSF phase is nothing but a limiting case of the VSF pre¬ 
viously described in the ladder case, when the vortices are 
of the size of a single plaquette in the square lattice. In 
addition, the CSF has nonvanishing density modulations. 
Were the Peierls phase density-independent, all the su¬ 
perfluid region would have nonvanishing chiral order, and 
the density modulations would disappear (not shown). 
Thus, the density-dependence in the Peierls phase has 
the effect of inducing finite density modulations and re¬ 
ducing the region with nonvanishing chiral order to that 
of p > 1. 


The phase transitions are in all cases found to be of 
second order, signaled by discontinuities in the second 
order derivative of the energy with respect to the chem¬ 
ical potential. The continuous vanishing of the BEC or¬ 
der parameter (SF-MI transition) and the bond-chiral 
order 77 (CSF-SF transition) also supports this assump¬ 
tion. Moreover, the phase diagram remains stable un¬ 
der increasing of the cluster size, as the CBMFT-4x2 
includes minor quantitative corrections to the phase bor¬ 
ders of the CBMFT-2x2. In particular, the CSF-SF 
phase boundary obtained with 2x2 and 2x4 basically 
overlap (Figs. [5] (a) and (b)). 


Comparing the phase diagrams for the ladder (Fig. [2| 
and the 2D square lattice (Fig. [5|, we observe that the 
modulated SF can be considered as the bulk counterpart 
of the MSF appearing in the ladder geometry. 


V. DYNAMICALLY PROBING THE 
DENSITY-DEPENDENT FIELD 


DDSM results in an intriguing dynamics that may be 
easily probed experimentally. We illustrate this point 
with the particular case of the ladder model (12) with 
0 1 = —02 = 0 and t x = t y . We are interested in the dy¬ 
namics of a defect (formed by either a doubly-occupied 
site, i.e. a doublon, or an empty site, i.e. a holon) cre¬ 
ated in a MI with p— 1, initially at site (1, j = 0). Note 
that this initial condition is chosen for simplicity of the 
analysis. The initial doublon or holon may be created 
in a more delocalized region of the ladder. The relevant 
conclusions about the expansion dynamics would be un¬ 
affected. Similar dynamics has been studied recently in 
the context of Bose Hubbard models without gauge fields 
[58l and may be observed in experiments with single site 
resolution 


For U t x quantum (particle/hole) fluctuations of 
the MI are irrelevant, and the defect expansion is like 
that of a single particle with a hopping t x (2 t x e ±1( ^) for 
the holon (doublon). Both holon and doublon propagate 
ballistically along the ladder, i.e. A j(r) = yj(j 2 )(r) ~ 
qr (we consider below the time r in units of h/t x for 
holons and h/2t x for doublons). The expansion coeffi¬ 
cient 7 is however markedly different. Holons do not 
experience any magnetic flux, and thus they propagate 
with a 0-independent 7 — \[2. In contrast, doublons 
experience a flux 20 and their trajectories are partially 
diverted by cyclotron motion. Hence 7 decreases with 
0 (Fig. [6]). The inset of Fig. [6] depicts examples of A j{r) 
for different 0. This situation has to be contrasted with 
the case of density independent magnetic fields. Here 
holons and doublons will both experience the same mag¬ 
netic flux 20 and propagate - up to a factor 2 due to 
bosonic enhancement - in the same way. 
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For lower U/t quantum fluctuations become relevant 
altering the defect expansion in an intriguing way. A 
perturbative treatment of the role of particle-hole fluc¬ 
tuations offers an instructive starting point of the study. 
Up to second-order one virtual doublon-holon pair may 
be created and annihilated which mediate new hoppings 
of the initial holon (doublon) of the form: 


_O f 

%(2) = -7^ENi + 2 ,i)0'.i|+H.c.] 


u 

2 txty 

u 


^2\Pi\j + l,k^i)(j,i\+U.c.], (21) 




where | i,j) denotes a defect at site (z, j), ai = 
e i 0 i p. = i -|_ e '<t>i (i 4 , for dou- 

blons (holons). In order to study the influence of quan¬ 
tum fluctuations beyond perturbation theory we perform 
t-DMRG calculations 03, with system sizes up to 100 
rungs keeping up to 1000 matrix states. 

As shown in Fig. [ 6 ] for U/t x = 10, for = 0 fluctua¬ 
tions speed up defect expansion; the expansion coefficient 
may reach values 7 ~ 1.6. This is intuitively clear since 
there are more processes expanding the defect along the 
ladder. This remains true for small 0. However, the pe¬ 


culiar phase dependence of the extra terms ( 21 ), modifies 


as well the effective magnetic flux experienced by the dou- 
blons. Indeed, for a sufficient large 0, fluctuations slow 
down the doublon expansion, i.e. they strengthen the 
cyclotron motion diverting the doublon expansion, cor¬ 
responding to an increase of the effective magnetic field 
experienced by the doublons. Moreover, quantum fluctu¬ 



ations make holon expansion ^-dependent due to virtual 
doublons. For sufficiently large </>, fluctuations slow down 
the holon expansion, i.e. holons experience an effective 
cyclotron motion induced by quantum fluctuations of the 
MI substrate. 


VI. SUMMARY 


Raman-assisted hopping may be used to induce 
density-dependent synthetic magnetism in cold lattice 
gases. In one dimensional systems this results in the 
interesting possibility of studying the anyon model [35 . . 
For ladders and 2D square lattices we have shown that 
these fields lead to a rich ground-state physics character¬ 
ized by the non-trivial interplay between density modula¬ 
tions and chirality. In two-leg ladders it is characterized 
by a density-driven Meissner- to vortex-superfluid transi¬ 
tion. Moreover, DDSM significantly affects the dynamics 
of particles in the lattice, leading in particular to an in¬ 
triguing expansion dynamics for doublons and holons in 
a MI, which presents a remarkable dependence on quan¬ 
tum fluctuations and may be used to reveal experimen¬ 
tally the DDSM. 

Although we have focused on ladders and 2 D square 
lattices, similar ideas may be applied to more general 
lattices, opening interesting possibilities for the realiza¬ 
tion of density-induced geometric frustration. In this 
work we discussed bosonic particles in the presence 
of DDSM. In Ref. [35 it is shown that also fermionic 
species may be a useful candidate for the realization of 
DDSM in cold atom experiments, since here the spurious 
(vi)-(ix) processes of section II E identically vanish. 
While in one dimensional systems this can be exploited 
to study the anyon Hubbard model, in two and higher 
dimension a significantly different model is realized. 
These possibilities will be examined in forthcoming 
works. 
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FIG. 6 : Expansion coefficient 7 of a defect along the lad¬ 
der as a function of the effective flux </> for t x — t y , and 
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symbols denote the t-DMRG results for the holon (doublon) 
expansion. Dashed (solid) curves denote single-particle (exact 
diagonalization) results for holons (doublons), which match 
well with the t-DMRG results for large U/t x . The inset 
depicts typical linear expansions of A j(r) for a doublon at 
U 00 and </>/7 T = 0 (solid line), 0.5 (dashed line) and 1 
(dotted line). 
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